Pakker

library(tidyverse)
library(janitor)
library(sf)
library(tmap)

Det er nemt at visualisere spatialt data

… medmindre man ikke aner, hvordan man gør!

Din nye bedste ven: sf

  • Den typiske måde at lege med spatialt data på er vha. GIS (Geographic Information Systems)-værktøjer designet til det

    • QGis, ArcGIS osv.

  • Programmer som R er dog løbende blevet udvidet med pakker, der gør det muligt at klare alting i det samme stykke software, som man bruger til andre ting

  • Det er dobbelt smart, fordi man har alting ét sted og bygget op omkring kode, der kan ændres og opdateres

  • Med sf er det blevet smooth sailing. Den måde, pakken håndterer det spatiale aspekt af et datasæt gør, at det ligner alle andre datasæt til forveksling

  • I kan kende næsten alle sf-funktioner på deres st_*-prefix (fx st_cast())

Verden før sf

Data ifølge {sp}

Data ifølge {sp}

Læs og gem data

Spatialt data opbevares typisk enten som

  1. Rå-data (ikke-spatialt format) men med x, y-koordinater
  • For alt andet end punktdata, kan det godt være lidt irriterende at arbejde med, fordi det så skal “samles” korrekt. Sidstnævnte støder I heldigvis sjældent på
  1. geojson
  • Har jeg faktisk aldrig brugt til lokale filer, kun til at læse ind online
  1. shapefile
  • Faktisk et lidt misvisende navn, fordi shapefiles består af flere filer
  • Jeg vil anbefale at opbevare hver “fil” (sæt af filer) i en separat mappe med samme navn som filerne

I kan læse geojson-filer med `sf::read_sf``:

url <- 
  "https://api.dataforsyningen.dk/kommuner?format=geojson"

kommuner_raw <- 
  read_sf(url)

rm(kommuner_raw)

Til shapefiles foretrækker (plejer?) jeg at bruge sf::st_read(). Her angiver dsn den mappe, jeres (flere!) filer ligger i. layer angiver navnet, de forskellige filer har til fælles (dvs. alt undtagen endelser, fx .shp)

europe <- st_read(dsn = "data/europe",
                  layer = "europe")
## Reading layer `europe' from data source 
##   `C:\Users\jevi\OneDrive - Epinion\Documents\GitHub\crashcourse_geoviz\data\europe' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 48 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -10.39023 ymin: 36.02593 xmax: 41 ymax: 71.17769
## Geodetic CRS:  WGS 84

I kan gemme jeres shapefiles ved hjælp af sf::st_write():

st_write(europe,
         dsn = "data/europe",
         layer = "europe",
         encoding = "UTF-8",
         delete_layer = TRUE,
         factorsAsCharacter = TRUE,
         driver = "ESRI Shapefile")
## Deleting layer `europe' using driver `ESRI Shapefile'
## Writing layer `europe' to data source `data/europe' using driver `ESRI Shapefile'
## Writing 48 features with 3 fields and geometry type Multi Polygon.

Where the magic happens: df$geometry

Det smarte ved sf’s måde at opbevare spatialt data på er, at ALT det spatiale ligger i den liste, der hedder geometry. Den er også “fredet”, og man kan ikke slette den med select(). Hvis man vil af med den (og det vil man af og til, det vender vi tilbage til!), er der en specifik måde til det.

Udover geometry ligner spatiale datasæt en tibble som I kender dem. Bemærk, at sf tilføjer geografisk metadata, når I printer datasættet (inkl. CRS, datatype og geografisk omfang, bounding box)

europe
## Simple feature collection with 48 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -10.39023 ymin: 36.02593 xmax: 41 ymax: 71.17769
## Geodetic CRS:  WGS 84
## First 10 features:
##             cntry           unit          sub_re                       geometry
## 1         Vatican        Vatican Southern Europe MULTIPOLYGON (((12.43838 41...
## 2  United Kingdom         Jersey Northern Europe MULTIPOLYGON (((-2.082227 4...
## 3  United Kingdom       Guernsey Northern Europe MULTIPOLYGON (((-2.520898 4...
## 4  United Kingdom    Isle of Man Northern Europe MULTIPOLYGON (((-4.392285 5...
## 5  United Kingdom United Kingdom Northern Europe MULTIPOLYGON (((-2.539355 5...
## 6         Ukraine        Ukraine  Eastern Europe MULTIPOLYGON (((38.20586 47...
## 7     Switzerland    Switzerland  Western Europe MULTIPOLYGON (((9.35 47.598...
## 8          Sweden         Sweden Northern Europe MULTIPOLYGON (((18.95645 57...
## 9           Spain          Spain Southern Europe MULTIPOLYGON (((1.592676 38...
## 10       Slovakia       Slovakia  Eastern Europe MULTIPOLYGON (((22.47305 49...

Der er intet mystisk eller unikt ved “spatialt” data

tribble(~name,~n_champs,~lat, ~lon,
        "Vejle Boldklub",5,55.71372762470475, 9.556222451831141,
        "Aarhus Gymnastikforening",5,56.13193948661876, 10.196519237914524,
        "AC Horsens",0,55.87162048058711, 9.857707600274574)
## # A tibble: 3 x 4
##   name                     n_champs   lat   lon
##   <chr>                       <dbl> <dbl> <dbl>
## 1 Vejle Boldklub                  5  55.7  9.56
## 2 Aarhus Gymnastikforening        5  56.1 10.2 
## 3 AC Horsens                      0  55.9  9.86
df <- tribble(~name,~n_champs,~lat, ~lon,
              "Vejle Boldklub",5,55.71372762470475, 9.556222451831141,
              "Aarhus Gymnastikforening",5,56.13193948661876, 10.196519237914524,
              "AC Horsens",0,55.87162048058711, 9.857707600274574)

df <- df %>% 
  st_as_sf(coords = c("lon", "lat"),
           crs = 4326)

tmap_mode("view")

tm_shape(df) +
  tm_bubbles(size = 2,
             col = "n_champs")

Datatyper

points

urbana <- st_read(dsn = "data/urbana",
                  layer = "urbana")
class(urbana$geometry)
## [1] "sfc_POINT" "sfc"
ggplot() +
  geom_sf(data = urbana,
          size = .01) +
  theme_void()

lines

rivers <- st_read(dsn = "data/rivers",
                  layer = "rivers")
class(rivers$geometry)
## [1] "sfc_MULTILINESTRING" "sfc"
ggplot() +
  geom_sf(data = rivers) +
  theme_void()

polygons

europe <- st_read(dsn = "data/europe",
                  layer = "europe")
class(europe$geometry)
## [1] "sfc_MULTIPOLYGON" "sfc"
ggplot() +
  geom_sf(data = europe) +
  theme_void()

Visualisér med ggplot2

Vi vil gerne vise størrelsen af urbane centre i Europa

Vi har også en idé om, at mange større byer ligger tæt på floder, fordi det historisk set har været en fordel at være tæt på vand, fordi det gjorde handel lettere.

Lad os starte med at vise kortet over Europa!

Fordi vi bruger sf-pakken til det hele, kan vi bruge tidyverse-funktioner til at manipulere data med – og ggplot2 til at visualisere det!

Vi har ikke brug for at vise koordinater hen ad akserne, så vi bruger theme_void():

ggplot() +
  geom_sf(data = europe) +
  theme_void()

Ikke dumt! Lad os tilføje vores data over byer. Det gør vi let ved at tilføje endnu en “geometry” til vores `ggplot:

ggplot() +
  geom_sf(data = europe) +
  geom_sf(data = urbana) +
  theme_void()

Aha! Det er jo meget fedt, at vi har data på byer rundt i hele verden, men lige her gør det nok mere skade end gavn.

Husk, at tidyverse-funktioner kan bruges på sf-datasæt, akkurat som på klassiske datasæt. Hvis der var en regionsvariabel i bydatasættet, kunne vi derfor let filtrere alt uden for Europa væk vha. dplyr::filter(). Det er der desværre ikke, så vi må ty til andre tricks.

Heldigvis har sf-pakken en lang række værktøjer, der hjælper os med at manipulere data. Fx kan vi bruge st_intersection() til at beskære ét datasæt med et andet:

eu_cities <- urbana %>% 
  st_intersection(.,
                  europe)

ggplot() +
  geom_sf(data = europe) +
  geom_sf(data = eu_cities,
          size = .8) +
  theme_void()

Meget bedre! Men, det bliver en smule tætpakket med alle de prikker. Lad os lege lidt med det visuelle:

ggplot() +
  geom_sf(data = europe, fill = "white") +
  geom_sf(data = eu_cities,
          size = 1,
          alpha = .5) +
  theme_void()

Da vi beskar bydatasættet med shapefilen med Europa-kortet, tilføjede den automatisk data fra sidstnævnte.

eu_cities
## Simple feature collection with 3424 features and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -9.251242 ymin: 36.1835 xmax: 40.98619 ymax: 69.74173
## Geodetic CRS:  WGS 84
## First 10 features:
##      scalerank featurecla area_sqkm min_zoom  tmp          cntry           unit
## 3294         9 Urban area     8.936      7.7 3294 United Kingdom         Jersey
## 3293         8 Urban area    34.616      7.6 3293 United Kingdom       Guernsey
## 3231         7 Urban area    18.935      7.0 3231 United Kingdom    Isle of Man
## 3198         9 Urban area    10.900      7.7 3198 United Kingdom United Kingdom
## 3199         9 Urban area    11.659      7.7 3199 United Kingdom United Kingdom
## 3201         9 Urban area    10.329      7.7 3201 United Kingdom United Kingdom
## 3202         7 Urban area    80.259      7.0 3202 United Kingdom United Kingdom
## 3203         9 Urban area     6.911      7.7 3203 United Kingdom United Kingdom
## 3204         8 Urban area    32.645      7.6 3204 United Kingdom United Kingdom
## 3205         9 Urban area    10.634      7.7 3205 United Kingdom United Kingdom
##               sub_re                   geometry
## 3294 Northern Europe POINT (-2.098328 49.19537)
## 3293 Northern Europe POINT (-2.550666 49.47685)
## 3231 Northern Europe  POINT (-4.48212 54.16456)
## 3198 Northern Europe POINT (-3.312512 57.65061)
## 3199 Northern Europe POINT (-3.635281 57.62275)
## 3201 Northern Europe POINT (-4.213206 57.47909)
## 3202 Northern Europe  POINT (-2.128716 57.1617)
## 3203 Northern Europe POINT (-2.467238 56.73593)
## 3204 Northern Europe POINT (-2.607162 56.59364)
## 3205 Northern Europe POINT (-2.738971 56.52185)

Det kan vi bruge til at tilføje flere farver til vores kort!

ggplot() +
  geom_sf(data = europe, fill = "white") +
  geom_sf(data = eu_cities,
          aes(fill = sub_re),
          shape = 21,
          size = 1,
          alpha = .5) +
  theme_void()

Ikke dumt! Men det er stadig svært at udlede så meget substantielt. Lad os prøve noget nyt: bruge de data, vi fik smidt på byerne før til at summere nogle mere aggregerede mål. På den måde kan vi visualisere landene (polygoner) med information fra byerne (punkter).

Bydatasættet har ikke data på population, men der er et mål for areal i kvadratmeter. Lad os antage, at det er en fin proxy for befolkning (her antager vi konstant befolkningstæthed, hvilket selvfølgelig ikke er korrekt, men lad os køre med det).

Her har vi to spatiale datasæt:

  1. Et over lande (polygoner)
  2. Et over byer (punkter), med lande- og regionsvariable fra det første

Det kan selvfølgelig kombineres! Igen, vi kan bruge alle de tidyverse-tricks vi kender og elsker.

Fordi vi i sidste ende gerne vil visualisere polygonerne, har vi ikke brug for det spatiale element i bydatasættet. Det skal bare merges igen. Derfor sletter vi geometry-listen. Hvorfor ?

  1. sf vil ikke merge to spatiale datasæt med en klassisk *_join(). Den vil meget hellere joine på det spatiale, men det har vi ikke brug for her
  2. Spatialt data er tungt. Som tommelfingerregel bør I derfor altid droppe geometry’en, når I ikke skal bruge den. Lad mig vise det nedenfor:
library(tictoc)

# aggreger vores areal-variabel inden for lande med spatialt data
tic()
city_cntry <- eu_cities %>% 
  group_by(cntry) %>% 
  summarize(area = sum(area_sqkm)) %>% 
  ungroup()
toc()
## 1.8 sec elapsed
# aggreger vores areal-variabel inden for lande med ikke-spatialt data
temp <- eu_cities %>% 
  st_drop_geometry()

tic()
city_cntry <- temp %>% 
  group_by(cntry) %>% 
  summarize(area = sum(area_sqkm)) %>% 
  ungroup()
toc() ; rm(temp)
## 0 sec elapsed

Forskellen er jo ikke enorm, men den er klar! Hvis I arbejder med større datasæt, kan det hurtigt blive en enorm fordel at tænke over, hvornår det spatiale aspekt af data mest bare er en hæmsko.

Nå, I digress. Tilbage til visualiseringen. Det nye city_cntry data kan vi let joine med vores landepolygoner med en simpel left_join(), og så er vi back in business!

df_europe <- europe %>% 
  group_by(cntry) %>% 
  summarize()

df_europe <- df_europe %>% 
  tidylog::left_join(.,
                     city_cntry,
                     by = "cntry")

Det gik jo smooth! Vi kan dog se, at tre lande (i polygondatasættet) ikke blev merged med noget data. Det er sådan noget man altid lige skal tjekke (og derfor jeg ELSKER {tidylog}-pakken), men her er det bare fordi der ikke var nogle større byer. Lad os for god ordens skyld bruge den information til at erstatte NA med 0 i vores nye variabel:

df_europe <- df_europe %>% 
  mutate(area = ifelse(is.na(area),
                       0,
                       area))

Nu kan vi lave vores nye kort!

ggplot() +
  geom_sf(data = df_europe,
          aes(fill = area)) +
  theme_void()

Bum! Det var lige det, vi gerne ville have. Men faktisk er det lidt dumt, at Rusland er med på kortet, når det ikke rigtigt er en del af Europa – og fordi vi alligevel kun viser en lille del af det. Lad os fjerne det med en simpel filter():

df_europe <- df_europe %>% 
  filter(cntry != "Russia")

ggplot() +
  geom_sf(data = df_europe,
          aes(fill = area)) +
  theme_void()

Sådan! Nu vi er i gang, kan vi godt prøve at gøre det lidt pænt.

ggplot() +
  geom_sf(data = df_europe,
          aes(fill = area),
          color = "white",
          size = .01) +
  scale_fill_viridis_c(name = "Sum of Area (km^2)\n ofUrban Areas") +
  theme_void() +
  guides(fill = guide_colorbar(direction = "vertical",
                               barheight = 12,
                               barwidth = .5))

Det eneste, der måske kan være lidt misvisende er, at vores befolkningsmål ikke er korrigeret for landenes areal. Det er heldigvis nemt at implementere med st_area() fra sf-pakken og lidt datarwrangling:

df_europe <- df_europe %>% 
  mutate(poly_area = st_area(.)) %>% 
  mutate(poly_area = unclass(poly_area))

df_europe <- df_europe %>% 
  mutate(relative_area = area/poly_area)

ggplot() +
  geom_sf(data = df_europe,
          aes(fill = relative_area),
          color = "white",
          size = .01) +
  scale_fill_viridis_c(name = "Urban Area as Share\n of Total Area") +
  theme_void() +
  guides(fill = guide_colorbar(direction = "vertical",
                               barheight = 12,
                               barwidth = .5))

Det tegner jo et lidt andet mønster.

Lad os, for at vende tilbage til den oprindelige analyse, se om datasættet over floder kan bidrage lidt!

# beskær data til kun at dække Europa
rivers <- rivers %>% 
  st_intersection(.,
                  df_europe)


ggplot() +
  geom_sf(data = df_europe,
          aes(fill = relative_area),
          color = "white",
          size = .01) +
  scale_fill_viridis_c(name = "Urban Area as Share\n of Total Area") +
  geom_sf(data = rivers,
          color = "orange") +
  theme_void() +
  guides(fill = guide_colorbar(direction = "vertical",
                               barheight = 12,
                               barwidth = .5))

ggsave(plot = last_plot(),
       filename = "plots/europe_static.png")

Visualisér med tmap

Til de fleste formål vil ggplot2 være bedre til at visualisere jeres spatiale data på. Syntaxen er også noget mindre logisk, selvom den overordnede struktur minder om ggplot2’s.

tmap har dog en central fordel: det er eminent til at lave interaktive kort!

Af samme grund bruger jeg det flittigt, når jeg arbejder med/manipulerer data, fordi det gør det nemt at inspicere datakilder og -kvalitet:

# sæt til interaktiv visning
tmap_mode("view")

tm_shape(df_europe) +
  tm_polygons(col = "relative_area", alpha = 0.85,
              border.col = "white", lwd = .8) +
  tm_shape(rivers) +
  tm_lines(col = "steelblue")

Bemærk i øvrigt, at i modsætning til geom_sf() i ggplot2, skal tmap have specificeret typen af spatial data (med andre ord, hvilken geometry, der er tale om).

html <- tm_shape(df_europe) + 
  tm_polygons(col = "relative_area", alpha = 0.85,
              border.col = "white", lwd = .8) +
  tm_shape(rivers) +
  tm_lines(col = "steelblue")

tmap::tmap_save(html,
                "plots/europe_interactive.html")

CRS: Projektioner af kloden

“All models are wrong, but some are useful”
– George Box

Vi vil ofte gerne vise en globe i et 2D-rum. Problemet er, at er fysisk umuligt. Det kan simpelthen ikke lade sig gøre. Det betyder, at alle kort I nogensinde har set (glober undtaget) er forkerte på én eller anden måde. Tror I mig ikke? Se her!.

Det er her, CRS’er (Coordinate Reference Systems) kommer ind i billedet. For at “løse” det helt fundamentale problem (og fordi glober er ret besværlige at gå rundt med) er der udviklet en lang række projektioner. Vi bruger med andre ord disse matematiske projektioner til at projicere en globe (oftest vores egen klode) i et 2D-rum på en måde, der fungerer til formålet.

Alle disse projektioner er “forkerte” på en eller anden måde, og der er altid en eller anden forvrængning af størrelse, form eller retninger. Men, de er jo ikke udviklet for sjov. De er udviklet, fordi vi har brug for dem, og fordi de er brugbare. De er ikke allesammen lige brugbare, og de er sjældent lige brugbare til samme formål. Med en let omskrivning af ovenstående citat kan vi beskrive dem som:

→ All projections are wrong, but some are useful

Det nok mest berømte eksempel på en projektion (og på, at de aldrig er perfekte…) er Mercator-projektionen (EPSG 4326 / 3857, se mere her). Den er berømt og berygtet for at repræsentere visse dele af kloden (især Afrika og Sydamerika) som værende mindre, end de faktisk er. Det er jo uden tvivl rigtigt nok. Se billedet nedenfor, hvor cirklerne angiver areal relativt til faktisk areal. Den mere ukvalificerede del af kritikken går så langt som til at fremhæve, at den er “forkert”. Det er den, men det er alle alternativer også.

Der er også spekuleret meget i, at projektionen skulle være udviklet med henblik på mentalt at (re)producere europæisk dominans, fordi den fremstiller Afrika som langt mindre relativt til Europa end det faktisk er:

“These maps of Africa, drawn up by a small group of western cartographers, symbolically reinforced Europeans’ sense of control over their mapped territories and subjects, but they didn’t betray much in the way of real information. Though they would have been seen as objective and impartial at the time, in retrospect it is clear how subjective, ideologically driven, and, in many ways, fantastical they were.”
– James Wan

Det er noget sludder. Det er også et ulogisk argument. I det 15. århundrede, hvor Mercator-projektionen blev udviklet, havde europæerne ikke behov for at fremstille det afrikanske argument som værende “småt” for at legitimere kolonialisering. Der var masser af “legitimitet” at hente qua religion, “civilisation”, vanvittige idéer om “raceforskelle” osv. Tværtimod havde europæiske magthavere, om noget, interesse i at fremstille deres oversøiske besiddelser som værende så store og imponerende som muligt.

Mercator-projektionen er uden tvivl relateret til europæisk kolonialisme, men på en lidt mere indirekte måde. Den var nemlig helt eminent at navigere efter, især på meget lange maritime strækninger og med simpel teknologi. Som, fx, på turene over Atlanterhavet i det 15. århundrede. På disse længere ture giver Mercator-projektionen ikke helt den korteste rute, men det giver den mest bullet-proof kurs.

Hvilken CRS skal I bruge?

Det vigtigste er, at I bruger den samme CRS i alle jeres objekter. Et lille trick er at definere en værdi i starten, som I løbende kan referere til. På den måde skal I kun ændre den ét sted, hvis det bliver aktuelt

my_crs <- 4326

# df <- read_sf() %>% 
#   st_transform(4326)
# 
# map <- read_sf() %>% 
#   st_transform(3359)
# 
# st_crs(df) == st_crs(map)
# 
# df <- df %>%
#   st_transform(crs = my_crs)
# 
# map <- map %>%
#   st_transform(crs = my_crs)
# 
# st_crs(df) == st_crs(map)

For mere information om brug af CRS’er i R, se: